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Abstract 

An algorithm which either finds an nonzero integer vector m for given t real n-dimensional vectors Xj , • • • , x, such that 
x ( r m = or proves that no such integer vector with norm less than a given bound exists is presented in this paper. The 
cost of the algorithm is at most 0(n A + n 3 log A(X)) exact arithmetic operations in dimension n and the least Euclidean 
norm A(X) of such integer vectors. It matches the best complexity upper bound known for this problem. Experimental 
data show that the algorithm is better than an already existing algorithm in the literature. In application, the algorithm 
is used to get a complete method for finding the minimal polynomial of an unknown complex algebraic number from 
its approximation, which runs even faster than the corresponding Maple built-in function. 

Keywords: integer relation, PSLQ, HJLS, algebraic number, minimal polynomial 



1. Introduction 

Given a real vector x = (jci, • • ■ , x n ) T e M", say a nonzero vector m = (m\ , ■ • • , m„) T e Z" is an integer relation for 
x if x r m = 2" =I jc,7«, = 0. How to detect an integer relation for a given real vector is an old problem. This is solved, 
for instance, by the PSLQ algorithm |§1 that together with related lattice reduction schemes such as LLL fl4ll . was 
named one of ten "algorithms of the twentieth century" by the publication Computing in Science and Engineering (see 
0). This paper considers a generalization of the problem. Let xi, ■ • • , x, be t vectors in W, and denote (Xi, • • • , x,) by 
X. A simultaneous integer relation (SIR) for Xi , • ■ ■ , x, is a vector m e Z" \ {0} such that X T m = 0, i.e. x^m = for 
i = 1, ■ ■ • , t. For short, we also call m an SIR for X. An algorithm which either finds an SIR for t real n-dimensional 
vectors or proves that no SIR with norm less than a given bound exists is presented in this paper. 

When t = 1, the problem of detecting integer relations for one rational or real vector is quite old. For two numbers 
(a\,a2), the venerable Euclidean algorithm does the job by computing the ordinary continued fraction expansion of 
the real number 01/02- For n > 3, many detecting algorithms under the names generalized Euclidean algorithm 
and multidimensional continued fraction algorithm were proposed. We refer the reader to 111 lL [8[] for comprehensive 
surveys. Among these integer relation algorithms, the LLL-based HJLS algorithm and the PSLQ algorithm |8] 
have been used frequently. 



To authors' known, the first algorithm to detect SIRs for several real vectors (f > 2) was presented in 111 111 , in which 
J. Hastad, B. Just, J. C. Lagarias, and C. P. Schnorr not only presented the HJLS algorithm to find integer relations for 
one real vector, gave the first rigorous proof of a 'polynomial time' bound for a relation finding algorithm, but also 
proposed a simultaneous relations algorithm. Unfortunately HJLS has a serious drawback: it is extremely unstable 
numerically (see J7y8[]). 

In their draft llql . C. Rossner and C. P. Schnorr proposed an algorithm which computes for real vectors xi, X2 
simultaneous diophantine approximation to the plane spanned by the vectors xi, X2 by using a modified HJLS al- 
gorithm. It can be seen as a special case t = 2 of the aforementioned problem. But for the moment, it is still in a 
preliminary state with some open problems. 
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The PSLQ algorithm Jit] is now extensively used in Experimental Mathematics, with applications such as identifi- 
cation of multiple zeta constants, a new formula for n, quantum field theory and so on (see HI Sill])- PSLQ employs 
a numerically stable matrix reduction procedure, so it is numerically stable in contrast to other integer relation algo- 
rithms. Moreover, it can be generalized to the complex number field and the Hamiltonian quaternion number field, 
but the corresponding outputs are in Gaussian integer ring and Hamilton integer ring respectively. For example, PSLQ 
will output (I, I, -I) T for the complex vector (1 + /, 1 + 21, 2 + I) T , where / = V-T. The reason is that Hermite 
reduction in PSLQ produces some Gaussian integers in the reducing matrix (see Section l2~2l ). Thus PSLQ can not be 
used to detect SIRs (in Z") for several real vectors. 

An algorithm to detect SIRs for t real vectors is presented in this paper. It uses a technique, similar to that in HJLS, 
to construct the hyperplane matrix and a method, generalized from PSLQ, for matrix reduction. The algorithm either 
finds an SIR for X if one exists or proves that there are no SIRs for X of norm less than a given size. The cost of the 
algorithm is at most 0(n A + n 3 log A(X)) exact arithmetic operations to detect an SIR for X, where n is the dimension 
of the input real vectors and A(X) represents the least Euclidean norm of SIRs for X. Although the same theoretic 
complexity as obtained for the HJLS simultaneous relations algorithm is proved , experiments show that the algorithm 
in this paper often perfoms better in practice. Furthermore, in contrast to PSLQ, our algorithm can be applied to 
detect an integer relation in Z" (rather than in the Gaussian or Hamiltonian integer rings) for complex or Hamiltonian 
vectors. Consequentially, a complete method to find the minimal polynomial of an approximately complex algebraic 
number is obtained by applying our algorithm. 

Our main contributions in this paper are the following: 

• We present a new algorithm to detect SIRs for several real vectors and show that its complexity matches the 
best one known for this problem (HJLS simultaneous relation algorithm). 

• We implement our algorithm in Maple by two schemes. The one uses software floating arithmetic (multipreci- 
sion floating point arithmetic) in all steps, and the other partially uses software floating arithmetic and mainly 
uses hardware floating arithmetic. Then we report many experimental results, which shows that our algorithm 
is relevant. 

• We successfully apply our algorithm to find the minimal polynomial of an approximately complex algebraic 



number. This strategy is different from some known LLL-based methods, such as H13LI12H . and is for all complex 
algebraic numbers rather than mere for real agebraic numbers in 111511. We also present many experiments, 
which shows that this newly complete method is efficient and even better than the Maple built-in function 
PolynomialTools : -MinimalPolynomial. 

The reminder of this paper is organized as follows. In Section|2] both preliminaries and main results of this paper 
are presented. The cost of our algorithm is analyzed in Section[3] Some empirical studies, further discussions and an 
application of our algorithm are included in Section[4] 



2. The Main Algorithm 

2.1. Notations and Assumptions 

Throughout this paper, Z, R, and C stand for the sets of integers, real numbers, and complex numbers respectively. 
For c 6 R, denotes an arbitrary integer closest to c, i.e. |_cl = L c + hi- All vectors in this paper are column vectors, 
and will be denoted in bold. If x e R", then ||x||2 represents its Euclidean norm, i.e. ||x||2 = V(x, x), where (*, *) is 
the inner product of two vectors. We denote the n xn identity matrix by /„. Given a matrix A = (ay), we denote its 
transpose by A T , its trace by tr(A), its determinant by \A\, and its Frobenius norm by \\A\\f = (tr(AA r ))^ 2 = Q] a 2 p 1 ^ 2 . 
We say that a matrix A is lower trapezoidal if ay = for i < j. The group of n x n unimodular matrices with entries 
in Z are denoted by GL(n, Z). 
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In what follows we always suppose that Xi, • 
t < n. Obviously, every X; is nonzero. Let X e 



• ,Xj € R" are linearly independent, where x, = ■ • • ,x,.„) r and 
l nxt be the matrix (xi, ■ ■ • ,x,) and suppose that X satisfies 



X\jl-t+l x 2,n-l+l ' ' ' Xt,n-t+l 
Xl n-t+2 x 2,n-t+2 ' ' 1 x t,n-t+2 



x l.n 



X 2,l! 



±0. 



(1) 



If X G R" x ' does not satisfy (QJ we can always (since X has rank f) exchange some rows of X to produce X' - CX 
such that X' satisfying ([TJ, where C is an appropriate matrix in GL(n, Z). In this case, we detect an SIR for X' . If m 
is detected as an SIR for X', then C T m is an SIR for X. 



2.2. A Method to construct a Hyperplane Matrix 

Definition 2.1 (Hyperplane Matrix). Let X = (xi, • • • , X { ) s R" x ' . A hyperplane matrix with respect to X is any matrix 
H € R" x ("-') such that X T H = and the columns of H span X ± = {y e W : xj y = 0, i = 1 , • • • , ?}. 

Given X = (xi, ■ • • , x,) e R" x ' satisfying ([TJ. We now present a method to construct a hyperplane matrix for X. 
The basic idea is from HJLS iflUl . The same strategy was also used in PSLQ, based on a partial-sum-of-squares vector 
and a lower-quadrature matrix factorization, instead of Gram-Schmidt orthogonalization. 

Let bj , • • ■ , b„ form the standard basis of R", i.e. the z'-th entry of b,- is 1 and others are 0. Perform the process 
of standard Gram-Schmidt orthogonalization to Xi, • • • , X t , bi, • • • , b„ in turn producing Xj , • • ■ , x*, bf , ■ • • , b* . Note 
that, since X satisfies ([TJ, we have b* f+] = ■ ■ ■ = b* = 0. 

Define Hx to be the n X (« — t) matrix (bp • ■ • ,b*_ ( ). From the following lemma, Hx = (bj, • • • ,b*_,) is a 
hyperplane matrix with respect to X € R" x '. 

Lemma 2.2. Let X e R" xf and H x be as above. Then 

L H^Hx = /„-/. 

2. \\H X \\ F = ^h~=t. 

3. (xj, • • • , x*, Hx} is an orthogonal matrix. 

4. X T Hx —0, i.e. Hx is a hyperplane matrix ofX. 

5. Hx is a lower trapezoidal matrix and every diagonal element of Hx is nonzero. 

Proof. Since every two columns of Hx are orthogonal, part 1 follows. And part 2 follows from part 1. Let X* = 
(x*, ■ • • , x*) r . Obviously, (x*, • • • , x*, Hx) is an orthogonal matrix. From part 3 and standard Gram-Schmidt orthog- 
onalization we have X* T Hx = and X — X*Q respectively, where Q is an appropriate t x t invertible matrix. Thus 
X T Hx = Q T X* T Hx = and hence that part 4 follows. We now prove part 5. Denote the k-th element of b* by b*,. 

The diagonal elements of Hx are b* l for i — 1 , ■ • • ,n — t. Before normalizing b* we have b* u — 1 - Y?k=\ x *u ~ ^'j=i b*w 
and at the same time, 

1 i-l 

o#iib^ = (b;,K> = i-2<"-2^ 

Thus all the diagonal elements of Hx are nonzero. Now we only need to show that Hx is lower trapezoidal. From 
standard Gram-Schmidt orthogonalization, we can check that b* k = (b*,bp = holds for i > k. This completes the 
proof. □ 



1 Our assumption ( < n is based on the following fact: Any SIR m is in the orthogonal complement space of span(x[ , ■ ■ ■ , x ( ). Since Xi , ■ ■ ■ , x, 
are linearly independent vectors in K" we have ; < n. So if t = n, then the dimension of the linear space span(xi, ■ ■ ■ , x f ) is n, hence that there 
exists no simultaneous integer relations for xi , ■ ■ ■ , x,. 
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2.3. Generalized Hermite Reduction 

We now study how to reduce the hyperplane matrix Hx- First we recall (modified) Hermite reduction as presented 
infl. 

Definition 2.3 (Modified Hermite reduction). Let H = (hij) be a lower trapezoidal matrix with hjj + and set 
D :- I„. For i from 2 to n, and for / from i - 1 to 1 by step -1, set q := [hij/hjj]; then for k from 1 to n, replace d^ 
by d^ - qdj t k- We say DH is the modified Hermite reduction of H and D is the reducing matrix of H. 

If the entries of H are complex numbers, then q = [hij/hjj - ] may be a Gaussian integer. Thus for a complex vector, 
PSLQ can only gives a Gaussian integer relation. 

Hermite reduction is also presented in ||8|], and is equivalent to modified Hermite reduction for a lower triangular 
matrix H with hjj + 0. Both of the two equivalent reductions have the following properties: 

1. The reducing matrix D e GL(n, Z). 

2. For all k > i, the (modified) Hermite reduced matrix H' = (h' ( ) = DH satisfies \h' k ; | < \h' u \/2 = \hij\/2. 

Unfortunately, (modified) Hermite reduction is not suitable to detect SIRs any more because it does not deal with 
the last t - 2 rows of Hx when 2 < t < n. In order that the reduced and reducing matrices of Hx e R" x(n_,) satisfy the 
two properties above, we generalize the Hermite reduction as follows. 

Definition 2.4 (Generalized Hermite Reduction). Let H be a lower trapezoidal matrix with h } -j + and set D := /„, 
H' = (h'i j) := H. For i from 2 to «, for /' from min{/ - 1, n - t] by -1 to 1, q :- [h'.j/h'jj']; for k from 1 to /', 
h' ik = h' ik - qh'j k ; for k from 1 to n, d^ '■= - qdj^. For every two integers S\, Sz e {« — t + 1, • ■ ■ ,n} satisfying 
si < S2, h' „_, = and h' n _ t + 0, exchange the ij-th row and the .?2-th row of D. We call DH the generalized 
Hermite reduction of H and D the reducing matrix. 

Obviously, generalized Hermite reduction is equivalent to modified Hermite reduction when t — 1 . In addition, 
we can easily check that generalized Hermite reduction retains the two properties mentioned above when 1 < t < n. 

There are two main differences between the (modified) Hermite reduction and the generalized Hermite reduction. 
Firstly, the last t — 1 rows of H will also be reduced by the first n - t rows of H in the generalized Hermite reduction, 
while the (modified) Hermite reduction not. Secondly, generalized Hermite reduction exchanges the si-th row and the 
52-th row of D if s\ < S2, h Su „- t = and /i 52i „_, + hold. This implies that if /i„-/ + i,„-r = after generalized Hermite 
reduction then /z„- /+ 2,«-/ = • • • = h lh „- t = 0. 

2.4. The Algorithm Description 

Based on the method to construct the hyperplane matrix and the generalized Hermite reduction, an algorithm to 
detect SIR for real vectors is proposed as follows. 

Algorithm 1 (Simultaneous Integer Relation Detection). 

Input: (x i , • • ■ , x,) = X e R" x ' satisfying (Q]) and a parameter y > 2/ V3. 

Initiation. 

• Compute the hyperplane matrix Hx and set H := Hx, B := /„. 

• Reduce the hyperplane matrix H by the generalized Hermite reduction producing the reducing matrix D. Set 
X T := X T D \ H := DH, B := BD \ 

Iteration. 

1. Exchange. Let H = (hi J). Choose an integer r such that Y\hrj\ ^ "/\hij\ for 1 < i < n — t. Let 

a := h v , B :- h r + u , 

A:=h r+1 , r+u 5:= JpTW. ( j 

Define the permutation matrix R to be the identity matrix with the r and r + 1 rows exchanged. Update X T :- 
X T R, H :- RH, B :- BR. 
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columns r and r + 1 be 



Update H :=HQ. 



2. Corner. Let Q := /„_,. If r < n — t, then let the submatrix of Q consisting of the r-th and (r + l)-th rows of 

6/6 -A/6 
A/6 B/6 

3. Reduction. Reduce H by the generalized Hermite reduction producing D. Set X T := X T D l , H := DH, B :- 
BD l . 

4. Termination. Compute G := Then there exists no SIR whose Euclidean norm is less than G. Denote 
B = (Bi, • • ■ ,B„), where Bj e R". If X r B; = for some 1 < j < n, or fc„_t,«-( = then 

Output: the corresponding SIR for X. 

Remark 2.5. The description of Algorithm 1 is similar to PSLQ. But the biggest difference is that Algorithm 1 used 
the generalized Hermite reduction. The (modified) Hermite reduction used in PSLQ is not suitable to detect SIRs for 
several real vectors, as mentioned early. PSLQ may be viewed as a particular case of Algorithm 1 when t — 1 . 

Remark 2.6. Given a complex vector z = x + yl in C" where x, y e R" and I = V-T, finding an integer relation 
(in Z") for z is equivalent to finding an SIR for (x, y). Thus Algorithm 1 can be used. For instance, let z = (2 + 
37,4 + 91, 8 + 277, 16 + 817, 32 + 243/) r . For finding an integer relation for z, first let x = (2,4, 8, 16, 32) r and 
y = (3, 9, 27, 81, 243) r . Running Algorithm 1 with y = 1.16, and x, y as its input vectors gives an SIR (6, 7, -9, 2, 0) r 
for (x, y), which, of course, also is an integer relation for z. This is one of the biggest differences between Algorithm 
1 and PSLQ since PSLQ can only give a Gaussian integer relation in Z[7]" for z rather than an integer relation in Z". 

Remark 2.7. Generally, t real n-dimensional vectors may have 0, 1, or up to n — t linearly independent SIRs. One can 
follow the strategy in [8, Section 6] to find them. 

Theorem 2.8. Let X — (xi, ■ ■ • ,x r ) satisfy (fJ]) and AiX) be the least Euclidean norm of SIRs for X. Suppose there 
exists an SIR for X. Then 

1 . An SIR for X will appear as a column of B after no more than 



I -u 



\og(y n -A(X)) 



iterations in Algorithm 1. 

2. If after a number of iterations no SIR has yet appeared in a column of B, then there are no SIRs of norm less 
than the bound 1 . 

From this theorem, Algorithm 1 either finds an SIR m for given real vectors Xi, ■• ■ ,x, such that x^m = or 
proves that no small simultaneous integer relation with Euclidean norm less than 1/||//||f exists. 

Moreover, it can be proved that the norm of the SIR for X output by Algorithm 1 is no greater than y" ~'~ l A(X). 
This is an important property of Algorithm 1, and the proof is similar to that of Theorem 3 in |8|]. 

Corollary 2.9. IfX € R" xf has SIRs, then there exists a y such that Algorithm 1 can find an SIR for X in polynomial 
time 0(n 4 + n 3 log A(X)). 

Proof. Let y — 2. Then Algorithm 1 constructs an SIR for X in no more than 

(n - tf(n + t - 1) + (n - t)(n + t - 1) log A(X) 

iterations. Algorithm 1 takes 0(n - t) exact arithmetic operations per iteration, and hence that 0((n - t) 4 + (n — 
f) 3 log A(X)) exact arithmetic operations are enough to produce an SIR for X. Since t < n, the proof is complete. □ 

Remark 2.10. All conclusions above also hold for complex numbers with y > y2, but the outputs of the corresponding 
variation of Algorithm 1 are in Gaussian integer ring. 
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3. Proof of TheoremE^I 



Given X = (Xi, • • • ,X;) e R" x ', let Hx be the hyperplane matrix obtained by the method introduced in section |2~2 
and let Px = HxH^. By expanding this expression, it follows that Px — I n - Tli=\ x "j x * T ■ Let m e Z" be an SIR for X 
Then it can be seen that P^m = m and HPxIlf 
matrix Q, 

1 < UDmlb = \\DP x m\\ 2 < \\DP x \\ F \\m\\ 2 

= \\DHx\\ F \\m\\ 2 = \\DH x Q\\ F \\m\\ 2 



V« - t. For any matrix D e GL(n, Z) and (n - t)x(n - t) orthogonal 



(3) 



where \\DP x \\f = \\DH x \\f follows from P T X = P x and P\ = P X - From G}, the part 2 of Theoreml2T8lfollows. 
Let H(k) be the result after k iterations of Algorithm 1 . 

Lemma 3.1. Ifhjj(k) — Ofor some 1 < j < n — t and no smaller k, then j — n — t and an SIR for X must appear as a 
column of the matrix B. 

Proof. By the hypothesis on k, all diagonal elements of H(k - 1) are not zero. Now, suppose the r chosen in the 
Exchange step is not n - t. Since generalized Hermite reduction does not introduce any new zeros on the diagonal, 
and from the Exchange step and the Corner step, we have that no diagonal element of H{k) is zero. This contradicts 
the hypothesis on k, and hence that our assumption that r < n-t was false. Thus r - n-t after the (k - l)-th iteration. 

Next we show that there must be an SIR for X appeared as a column of the matrix B. We have X 7 Hx = from 
Lemma I2T21 and hence that = X T BB l H x = X T BB l H x Q = X T BH(k - 1), where Q is an appropriate orthogonal 



(n- t)x(n- 



t) matrix. Let (zj, 









X T B, where z,- — • • • ,Zi, n ) T - Then 



X T BH(k ■ 



1) 



Z"k=„-t Zl,khk,n-t(k ~ 1) 
£*=n-f Zt,khk,n-t(k ~ 1) 



H(k-l) 



Z\,n-th n -l.n-t(k ~ 1) 



Zt,n-th n -t,n-t(k — 1) 



(4) 



We know /i„_ (+ i.„_((A;- 1) = and h n - t n - t (k— 1) + from h n - ttn - t {K) = 0. From Definition l2~4l and h„-,+\ „_,(£- 1) 
we have h„-, + 2,„-t(k — 1) = • • • = h nM -,(k - 1) = which implies the last equality in (0). Since h n - t)n -t(k -1)^0, it 
follows that Zi, n -t - ' ' ' — Zt.n-t - 0. Thus the (« - f)-th column of B is an SIR for X. □ 



From the analysis above and Lemma BTl the correctness of Algorithm 1 is proved. From the iteration of Algorithm 
1, is decreasing with respect to k. Thus if there exist SIRs for X, Algorithm 1 can always find one. 

Definition 3.2 (II function). Let A(X) be the least norm of SIRs for X. For the £-th iteration in Algorithm 1, define 

Tl(k)= P[ min\y"-'A(X),- 

\<j<n-t 



l 



Lemma 3.3. For k > l we have 



1. (y"-'A(X)) 

2. n(k) > 



> Tl(k) > l. 



I). 



The routine of analyzing the number of iterations in JUl can be carried over here with redefining the II function as 
above, so we state Lemma [331 directly without proof. From this lemma, it follows that the II function is increasing 
with respect to k and has an upper bound for a fixed y e (2/ V3, +oo). From Definition [372] we can infer 11(0) > 1. 
And from Lemma [331 we know that 



(y"-'A(X)) 



> n(Ar) > 



4y 2 



\y 2 +4 



Solving k from this inequality gives the part 1 of Theorem l2.8l as was to be shown. 
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4. Empirical Study and Further Discussion 



4.1. Implementation 

All discussions above are based on exact arithmetic operation, i.e. uses the Blum-Shub-Smale model 01 of 
computation. The reason is that Algorithm 1 involves real numbers. Thus Algorithm 1 can only be implemented 
using floating point arithmetic on computer. Both Algorithm 1 and the HJLS simultaneous relations algorithm when 
t = 2, i.e. detecting an SIR for two real vectors, were implemented in Maple 13 under multiprecision floating 
point arithmetic (one level scheme). Like PSLQ (see [1, Section 5]), it is possible to perform most iterations using 
hardware floating point arithmetic, with only occasional depending on multiprecision arithmetic. So the two level 
implementation of Algorithm 1 is also developed by the first author. It partially uses software floating arithmetic and 
mainly uses hardware floating arithmetic (Cf. JH Section 5] for details). 

It is well known that the Gram-Schmidt orthogonalization algorithm is numerically unstable ficil . so in our im- 
plementations we construct the hyperplane matrix by using QR decomposition, instead of Gram-Schmidt orthogonal- 
ization. In contrast to HJLS, the iteration step in Algorithm 1 is not based on the Gram-Schmidt orthogonalization, 
but on LQ decomposition (this is equivalent to QR decomposition). Householder transformations are used in our 
implementations to compute these decompositions. Thus our implementations of Algorithm 1 is numerically stable. 



4.2. Experimental Result 

In theory, the cost of Algorithm 1 (in Corollary \2.9i matches the best complexity upper bound known for this 
problem (Cf. section 5]), whereas in practice Algorithm 1 usually needs fewer iterations. For Xi = (1 1, 27, 31) 
and X2 = (1,2, 3) r , HJLS outputs (19, -2, -5) T after 5 iterations while Algorithm 1 outputs the same SIR after only 3 
iterations. 



No. 


n 


itr H jLS 


itrsiRD 


tfiJLS 


tsiRD 


1 


4 


15 


8 


0.063 


0. 


2 


4 


13 


6 


0.062 


0. 


3 


4 


21 


11 


0.094 


0.015 


4 


5 


25 


12 


0.109 


0.016 


5 


5 


27 


7 


0.141 


0. 


6 


5 


21 


10 


0.094 


0. 


7 


30 


51 


7 


0.922 


0.125 


8 


54 


34 


9 


2.203 


0.453 


9 


79 


34 


5 


4.860 


0.625 


10 


97 


37 


5 


7.438 


1.047 



No. 


n 


itr H jLS 


itrsiRD 


tHJLS 


tsiRD 


11 


128 


45 


5 


13.765 


1.687 


12 


149 


29 


2 


19.016 


1.610 


13 


173 


26 


3 


26.812 


2.421 


14 


192 


29 


5 


34.218 


3.563 


15 


278 


28 


5 


85.797 


8.860 


16 


290 


35 


4 


95.656 


8.328 


17 


293 


23 


4 


98.062 


8.750 


18 


305 


22 


3 


109.187 


8.063 


19 


316 


19 


3 


120.187 


8.766 


20 


325 


18 


2 


129.031 


6.953 



Table 1 : Comparison of performance results for HJLS and Algorithm 1 

The purpose of the trials in TableQ]is to compare the performances of HJLS and Algorithm 1 when t = 2. All of 
the tests were run on AMD Athlon™ 7750 processor (2.70 GHz) with 2GB main memory. 

In Table [1] n gives the dimension of the relation vector, itrnjis and itrsiRD are the numbers of iterations of HJLS 
and Algorithm 1 respectively, and the columns headed tnjis and t$iRD give the CPU run time respectively of the two 
algorithms in seconds. The 20 trials in TableQ]were constructed by Maple's pseudo random number generator. The 
first 6 trials are for low dimension, and others for higher dimension. 

The results show that Algorithm 1 appears to be more effective than HJLS. In all 20 trials, the number of iterations 
of Algorithm 1 is less than that of HJLS. It is still true that Algorithm 1 usually needs fewer iterations than HJLS for 
more tests. This leads that the running time of Algorithm 1 is much less than HJLS. With the dimension n increasing, 
the difference between the efficiencies of Algorithm 1 and HJLS is increasingly notable. On average, the running time 
of Algorithm 1 is less than 1/10 (based on the data in Table Q]) of the running time of HJLS. 



: The package is available from http : // Cid-5dbbl6a211c63a9b . skydrive . live . com/self . aspx/ .Public/sird. rar 
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To some extent, the number of iterations and the cost of Algorithm 1 are related to the parameter y. In practice we 
have found that for many examples, larger values of y are more effective in finding SIRs for X. For Xi = (86, 6, 8, 673) r 
and X2 = (83, 5, 87, 91) r , if we choose y — 2 then Algorithm 1 outputs (-32, -747, 63, 10) r after 10 iterations, how- 
ever, if we choose y = 93, Algorithm 1 outputs (-35, -2624, 157, 26) r after only 7 iterations. It is worth mentioning 
that, in the example above, both of the two different output vectors are SIRs for (x;,X2) and they are linearly inde- 
pendent and hence that all SIRs for (xi,X2) can be obtained from them. Going on choosing a y larger than 93 in this 
example, after many tests, the authors find that the number of iterations is always 7, and it will not decrease any more. 
Based on this observation, all results in Table[T]are obtained under the condition that y = 1000. 

In general, a larger y requires a higher precision in authors' tests. Usually, a high precision leads a large height of 
the output SIR and a large cost of memory. If one sets y = 1.15470053838 (> 2/ V3J, for about 60% of our whole 
tests, the height of the vector returned by Algorithm 1 is less than that of HJLS simultaneous relations algorithm, 
however the number of iterations turns large. 

This means that we should try to find the balance between the number of iterations and the precision because both 
of them are relevant to the running time of Algorithm 1 . So in practice, what are the best choices for y needs further 
exploration. 

4.3. An Application 

We end this paper with an application of Algorithm 1 to find the minimal polynomial of a complex algebraic 
number from its approximation. 

Example 4.1. Let a — 2+ y/3I. We know that the minimal polynomial of a in Z[x] is 7-4X+X 2 . Let a = 2.000+1.732/ 
be an approximation to a with four significant digits. Let Vi = (1., 2., l.) r and V2 = (0., 1.732, 6.928) r be the real part 
and the imaginary part of (1, or, a 2 ) T respectively. Feeding Algorithm 1 with v\, \2 as its input vectors gives an SIR 
for Vi, V2 after 2 iterations. The corresponding matrices B are 



2 


1 







7 





2 


-1 


-1 







-4 





-1 








1 . 




1 


1 






It is obvious that the first column of the latter one is an SIR for Vi and V2 and corresponds to the coefficients of the 
minimal polynomial of a. However, if one takes only 3 significant digits for the same data, after 3 iterations Algorithm 
1 outputs (1213, -693, 173) r , which is an exact SIR for (l.,2., l.) r and (0., 1.73, 6.93) r , but does not correspond to 
the coefficients of the minimal polynomial of a. For this reason, how to appropriately control the error also is an 
interesting problem. 

Generally, for computing the minimal polynomial of an algebraic number a with degree n, we detect an integer 
relation for (1, a, ■ ■ ■ , a") T . If a e C, we detect an SIR for (1, Re(ar), • • ■ , Re(ff")) r and (0, Im(or), • ■ ■ , \m(a n )) T by 
Algorithm 1 under a proper decimal precision. The output vector corresponds to a polynomial of degree «, whose 
primitive part must be the minimal polynomial of a. 

As mentioned early, Algorithm 1 has been implemented in two schemes (one-level, two-level). Using our two level 
implementation of Algorithm 1, the authors obtain the following polynomial of degree 84 from an approximation to 
a = 3 1/6 - 2 1/7 7 with Digits: =1300 in Maple 13. It is easy to check that this polynomial is the exact minimal 
polynomial of a. 

5067001 + 783962907 x 36 + 21027764272536 x 40 - 7504504 x 42 
+ 836396 1 8394696 x 34 - 3668308 1 862336 x 38 + 770305668258672 x 32 
+ 142394998636968 x 28 + 1254656434122 x 30 + 1370000831472 x 10 
+ 34207465357611 x 12 - 284692059376032 x 14 + 24758 141678424 x 16 
+ 1909595 10258972 x 18 + 2306173886216928 x 20 + 99120704967648 x 22 
+ 641 11 149001809 x 24 - 3029254676588448 x 26 + 250312437648 x 44 
+ 2 1 89 1 87 x 48 - 1 1 26 1 5776 x 50 - 486486 x 54 + 8 1 08 1 x 60 - 378550368 x 2 
+ 1 1935794528 x 4 + 19043 11 10646 x 6 + 3293025660288 x 8 + x 84 
+ 88074554904 x 52 + 240 x 56 + 1041237288 x 58 + 1952496 x 64 - 9828 x 66 
+ 24 x 70 + 819 x 72 - 42 x 78 + 2212809521832 x 46 
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Table 2: Running times for the two implementations of Algorithm 1 



Some performance results are reported in Table [2] In this table, r and s (s is an odd integer number) define the 
constant a — 3^ r - 2 1/,J 7, which is an algebraic number of degree 2rs, and n - 2rs + 1. The column headed "Digits" 
gives a sufficient precision in decimal digits, while itr and t are the number of iterations and running times required 
for the correct output respectively, where the suffix SIRD is for one-level and TLSIRD for two-level. Every output 
vector in Table [^corresponds to the coefficients of the exact minimal polynomial of a's. 

From Table [2] it can be seen that the two-level program is much faster than the one-level program since the two- 
level program only involves multiprecision operation partially. The reason is that not all operations of Algorithm 1 
have to use high precision. As a matter of fact, the two-level program performs most of the iterations using IEEE 
hardware arithmetic. Thus the running times can be dramatically reduced. 

Using the case of t — 1 (it is PSLQ in fact) and t — 2 of Algorithm 1 for real and complex algebraic numbers, 
respectively, we get a new complete method (The method in 11 1 311 is only for real algebraic numbers.) to recover the 
minimal polynomial of an arbitrary algebraic number from its approximation and degree. It should be noted that since 
this method depends on integer relation detection that is based on a generalization of Euclidean algorithm [9], it is 
different from LLL-based algorithms, such as |Q~3, 121. 

In practice, the presented method is efficient. For a : = V21 + y43I, our procedure MiniPoly takes 1 .062 seconds 
for outputting the exact minimal polynomial of a, whereas the Maple built-in function MinimalPolynomial in 
PolynomialTools package that is LLL-based takes 6.032 seconds under the same decimal precision Digits : =500. 



5. Conclusion 

Using a method to construct a hyperplane matrix and the generalized Hermite reduction, a new SIRs detecting 
algorithm, Algorithm 1 , is presented in this paper. It runs faster than the HJLS simultaneous relation algorithm through 
the authors' Maple package. Applying the algorithm, we obtain a complete method to find the minimal polynomial 
of an approximately algebraic number, which is even faster than the corresponding Maple built-in function. 
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